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Fidelity susceptibility made simple: 

A unified quantum Monte Carlo approach 
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The fidelity susceptibility is a general purpose probe of phase transitions. With its origin in 
quantum information and in the differential geometry perspective of quantum states, the fidelity 
susceptibility can indicate the presence of a phase transition without prior knowledge of the local 
order parameter, as well as reveal the universal properties of a critical point. The wide applicability 
of the fidelity susceptibility to quantum many-body systems is, however, hindered by the limited 
computational tools to evaluate it. We present a generic, efficient, and elegant approach to compute 
the fidelity susceptibility of correlated fermions, bosons, and quantum spin systems in a broad 
range of quantum Monte Carlo methods. It can be applied both to the ground-state and non¬ 
zero temperature cases. The Monte Carlo estimator has a simple yet universal form, which can be 
efficiently evaluated in simulations. We demonstrate the power of this approach with applications 
to the Bose-Hubbard model, the spin-1/2 XXZ model, and use it to examine the hypothetical 
intermediate spin-liquid phase in the Hubbard model on the honeycomb lattice. 
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Phase transitions highlight the beauty of universal¬ 
ity, despite the great diversity of nature. For example, 
one finds a unified description for systems ranging from 
ultracold bosons [1, 2] to magnetic insulators [3-5] on 
the verge of a phase transition. Phase transitions origin 
from the competition between different tendencies when a 
macroscopic system tries to organize itself. Thermal fluc¬ 
tuations can drive classical phase transitions at non-zero 
temperatures, while quantum phase transitions can occur 
even at zero-temperature because of the competition be¬ 
tween non-commuting terms in the quantum mechanical 
Hamiltonian [6, 7]. At the phase transition point, phys¬ 
ical observables often exhibit singular behavior. In this 
respect, phase transitions are the most dramatic manifes¬ 
tation of the laws of statistical and quantum mechanics. 

Traditional descriptions of phase transitions are based 
on low-energy effective theories of local order parame¬ 
ters, which have had enormous success in explaining var¬ 
ious phase transitions of superfluids, superconductors [8], 
and quantum magnets [9, 10]. However, in recent years, 
exceptions to this Ginzburg-Landau-Wilson paradigm 
have emerged [11]. In particular, topological phase tran¬ 
sitions [12-14] do not have a local order parameters on 
either side of the phase transition. Therefore, new the¬ 
oretical tools are needed to search for and characterize 
these new quantum phases and phase transitions. Many 
concepts in quantum information science [15], such as 
the quantum fidelity and quantum entanglement, have 
proven to be useful [16, 17]. Having a point of view which 
is totally different from the traditional condensed matter 
approach, they do not assume the presence of a local 
order parameter and thus offer new perspectives of the 
phase transitions and their universalities. 

Specifically, we consider the following one-parameter 
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family of Hamiltonians with a driving parameter A, 

H{\) = H 0 + \Hi. (1) 

As A changes, the system may go through one or several 
phase transition(s) because of the competition between 
Ho and Hi. The quantum fidelity measures the distance 
on the manifold of A, which is defined as the overlap 
between the ground-state wavefunctions at two different 
values of the driving parameter, 


F(A,e)-=|(tf 0 (A)|¥o(A + e))|, (2) 


where H{ A) |T n (A)) = E n ( A) |T n (A)) and n = 0 corre¬ 
sponds to the ground state. Unless otherwise stated, we 
assume the wavefunctions are normalized and there is no 
ground-state degeneracy. It is anticipated that the fi¬ 
delity will exhibit a dip when the two wavefunctions are 
qualitatively different, e.g., when they belong to differ¬ 
ent phases [18]. This wavefunction overlap is also related 
to the Anderson’s orthogonality catastrophe [19] and the 
Loschmidt echo in quantum dynamics [20]. 

Since in general the quantum fidelity vanishes expo¬ 
nentially with the system size for a many-body system, 
it is more convenient to study the change of its logarithm 
with respect to the driving parameter, called the fidelity 
susceptibility [21], 


Xf( A) 


d 2 InF 


de 2 


e=0 


(3) 


The first-order derivative vanishes because F is at its 
maximum when e = 0. In general, the fidelity suscep¬ 
tibility is an extensive quantity away from the critical 
point, but it exhibits a maximum or even diverges at 
the critical point, thus indicating a quantum phase tran¬ 
sition [22, 23]. Similar to conventional thermodynamic 
quantities, it also follows a scaling law close to the criti¬ 
cal point [22-25], which can be used to extract universal 
information about the phase transition. An important 
feature of the fidelity susceptibility is that it can reveal 
a phase transition without prior knowledge of the local 
order parameter. This makes it suitable for detection 
of topological phase transitions [26-29] and for tackling 
challenging cases where an in-depth understanding of the 
underlying physics is still lacking [30, 31]. Interestingly, 
the fidelity susceptibility may also be accessible to exper¬ 
iments [32-34]. 

Despite its appealing features, the difficulty in cal¬ 
culating the fidelity susceptibility has hindered its use 
in numerical simulations. Many previous studies were 
thus limited to the cases where the ground-state wave- 
function overlap could be calculated from the analyt¬ 
ical solution, exact diagonalization, or density-matrix- 
renormalization-group (DMRG) methods [16]. 

There are several equivalent formulations of the fidelity 
susceptibility Eq. (3), which reveal different aspects of 


the quantity. From a computational point of view, they 
offer direct ways to calculate the fidelity susceptibility 
without the need to perform numerical derivatives of the 
fidelity as in Eq. (2). 

(a) Expanding |T 0 (A + e)) for small e, one can cast the 
definition Eq. (3) into an explicit form [22, 35], 


m _ <0A*O|0A*O> (^ol^A^o) <3 a*o|* 0> m 

XF[ ’ <*o|*o> <*o|*o> <*o|*o> ' 1 ’ 

The above form does not assume properly normalized 
wavefunctions |Tq). Equation (4) reveals the geometric 
content of the fidelity susceptibility [22, 35], since this 
expression is the real part of the quantum geometric ten¬ 
sor [36]. 

(b) Alternatively, one can calculate the first-order per¬ 
turbation for |Tq(A + e)) and get 


xf{ a) = 

n/0 


|(M/»(A)|g 1 |v|/o(A ))| 2 

K(A)-£ 0 (A)] 2 


(5) 


Compared to Eq. (4), Eq. (5) does not contain derivatives 
but involves all eigenstates and the full spectrum. It 
explicitly shows that xf > 0 and suggests the divergence 
of xf when the energy gap of the system closes. 

(c) Reference [21] views Eq. (5) as the zero-frequency 
component of a spectral representation, thus a Fourier 
transform is performed to obtain an alternative expres¬ 
sion, 


r c 

Xf( A) = / dr (^o\H 1 (r)H 1 \^o)-(^o\H 1 \^ o y 
Jo 


(6) 


where H\{r) = e Hr H\ c ~ Ht . Equation (6) has the form 
of a linear-response formula and is computationally more 
friendly than Eq. (4) or Eq. (5). References [24, 25] gen¬ 
eralize it to non-zero temperature by replacing the inte¬ 
gration limit with /3/2, 


Xf{ A) = 


rP /2 
/ dr 

Jo 


(H\ ( t ) Hi) — {HiY 


(7) 


where (...) denotes the thermal average at inverse tem¬ 
perature /3. Besides reducing to Eq. (6) as (3 —>> oo, 
Eq. (7) is nevertheless a well-defined quantity at nonzero 
temperatures. It bounds the divergence of an alterna¬ 
tive “mixed state” fidelity susceptibility [37, 38], which is 
based on the Uhlmann fidelity [39, 40], and both quanti¬ 
ties follow the same scaling law close to a quantum crit¬ 
ical point [24, 25]. In general, the evaluation of Eq. (7) 
is still a formidable computational task, which requires 
ad hoc implementation depending on the details of the 
Hamiltonian. For example, the fidelity susceptibility for 
two-dimensional quantum spin systems is calculated in 
Refs. [24, 25] using quantum Monte Carlo method; while 
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for a one-dimensional quantum spin system, Ref. [41] 
computed it using transfer-matrix DMRG method. The 
nontrivial implementation of these specific approaches 
and the overhead in the calculation still limits the wide 
applicability of the fidelity susceptibility approach to a 
broad range of quantum many-body systems. 

In this paper we present a simple yet generic approach 
to compute the fidelity susceptibility in a large vari¬ 
ety of modern quantum Monte Carlo methods, includ¬ 
ing the continuous-time worldline [42-45] and stochastic 
series expansion (SSE) [46] methods for bosons and quan¬ 
tum spins, and the diagrammatic determinantal methods 
for quantum impurity [47-50] and fermion lattice mod¬ 
els [51-53]. In all cases, the Monte Carlo estimator is 
generic and the implementations are straightforward. As 
long as the quantum Monte Carlo simulation is feasible 
(not hindered by the sign problem), the fidelity suscepti¬ 
bility can be easily calculated. Our finding can boost the 
investigation of quantum phase transitions from a quan¬ 
tum information perspective and becomes especially ad¬ 
vantageous for the exploration of exotic phases beyond 
the Ginzburg-Landau-Wilson paradigm. 

The organization of the paper is as follows. In Sec. II 
we present our estimator for the fidelity susceptibil¬ 
ity, and discuss its implementations in various quantum 
Monte Carlo methods. Section III presents derivations 
of the estimator. In Sec. IV we demonstrate the power 
of the fidelity susceptibility approach with applications 
to various models, including correlated bosons, fermions, 
and quantum spins, using a variety of quantum Monte 
Carlo methods. Section V discusses the relation between 
the zero-temperature and non-zero temperature estima¬ 
tors for the fidelity susceptibility and compares them to 
the previous approaches [24, 25]. We conclude with fu¬ 
ture prospects in Sec. VI. 


II. RESULTS 


We first present our results on the estimator of the fi¬ 
delity susceptibility in a general setting, then discuss its 
implementations in various QMC methods, including the 
continuous-time worldline [42-45] and diagrammatic de- 
terminantal approaches [47-53], and the stochastic series 
expansion method [46]. In all cases, the fidelity suscepti¬ 
bility can be measured with little effort. 


A. Universal Covariance Estimator 



Figure 1. Measurement of the fidelity susceptibility in the 
(a) non-zero temperature formalism and in the (b) ground- 
state projection scheme. Each red object represents a term 
in the driving Hamiltonian XH±, denoted as a vertex. To 
measure the fidelity susceptibility Eqs. (9,11), we divide the 
imaginary-time axis into two halves and count the number of 
vertices kL and kR respectively. The non-zero temperature 
formalism allows an arbitrary division because of the periodic 
boundary condition in the imaginary-time axis, while in the 
ground-state projection scheme the division has to be at /3/ 2. 


where the second summation runs over all the Monte 
Carlo configurations of a given expansion order k. The 
detailed meaning of the configuration depends on the spe¬ 
cific QMC algorithm and will be explained in the next 
subsection. Figure 1(a) depicts a generic configuration, 
where the k objects residing on the periodic imaginary- 
time axis represent the vertices XH\ in the expansion, 
with a Monte Carlo weight A k w(Ck) for this configura¬ 
tion. QMC simulations [42-53] sample the summation 
over k and C& on an equal footing. Specific algorithms 
differ by the detailed form of w(Ck) and by the sam¬ 
pling schemes. Nevertheless, these QMC methods share 
a unified framework provided by Eq. (8), which is the 
only requirement for the estimator of the fidelity suscep¬ 
tibility Eq. (7) to possess an appealing universal form in 
non-zero temperature QMC simulations, 


T/0 _ - (&l) ( kR ) . 

X F - 2A 2 , W 

where kL and kR are the number of vertices residing 
in the range [/3/2,/3) and [0,/3/2) of the imaginary-time 
axis, respectively, shown in Fig. 1(a). In practice, how¬ 
ever, because of the periodic boundary condition on the 
imaginary-time axis, the division of the time axis to 
halves may be done at an arbitrary location. Moreover, 
it is even possible to perform multiple measurements on 
the same configuration by generating several random di¬ 
visions. 

QMC methods [42-53] can also be utilized at zero 
temperature, where the unnormalized ground-state wave- 
function is obtained from an imaginary-time projection 


Many modern QMC methods [42-53] share a unified 
conceptual framework, namely that the partition func¬ 
tion is calculated as a perturbative series expansion, 

oo 

Z = Tr(e-e 6 ) ( 8 ) 

k =0 Ck 


|#o) = lim e"^ /2 I <S> T ). (10) 

/3 — yoo 

Here (3 is a projection parameter and the trial wavefunci- 
ton \S&t) shall not be orthogonal to the true ground state. 
A similar framework as Eq. (8) applies, except that one 
now samples from the overlap (Tt| e~^ H |4/t) instead of 
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in the form of Eq. (1) and perform a time-dependent ex¬ 
pansion in \H\, 


r P 
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dr k x 
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Figure 2. (a) Measurement of the fidelity susceptibility in 

a continuous-time worldline QMC simulation of bosons and 
quantum spins, where one counts the number kinks kL and kR 
after division of the imaginary-time axis, (b) In the diagram¬ 
matic determinantal QMC simulation of correlated fermions, 
the number of interaction vertices kL and kR are counted. 


the partition function. In the projection scheme, the fi¬ 
delity susceptibility has the estimator 


v t=o 
X F 


- (fc l) (kR) 

A 2 


(ii) 


where kL and kR are the number of vertices for the bra 
and ket states, which reside in the range [/3/2, /?) and 
[0, /3/2) of the imaginary-time axis respectively, shown in 
Fig. 1(b). Since the fidelity susceptibility is non-negative, 
the covariance formula Eq. (11) reveals positive correla¬ 
tion of kL and kR in a Monte Carlo simulation. 

Equation (9) and Eq. (11) are the central results of the 
paper. As is obvious from the discussions in this section, 
neither details of the Hamiltonian, nor the statistics of 
the system need to be specified. These estimators are 
thus general and can be readily implemented in a variety 
of QMC methods for correlated fermionic, bosonic, or 
quantum spin systems [42-53]. 


, ( 12 ) 


which obviously fits in the general framework of Eq. (8) . 
In the continuous-time worldline approach [42-45], the 
Hi term corresponds to hoppings of bosons or spin 
flips, depicted as kinks of the worldlines in Fig. 2(a). 
In continuous-time diagrammatic determinantal ap¬ 
proaches [47-53], XHi contains the fermion interactions, 
drawn as interaction vertices in Fig. 2(b). Equation (12) 
has the form of a grand canonical partition function for 
a classical gas, where A plays the role of fugacity and 
k is the number of certain classical objects (kinks or 
vertices) residing on the imaginary-time axis. Typical 
updates of continuous-time diagrammatic determinantal 
approaches [47-53] consists of randomly inserting or re¬ 
moving vertices, which are identical to the updates of 
grand canonical Monte Carlo method for molecular sim¬ 
ulations [54, 55]. For bosons and quantum spins there are 
more effective non-local updates such as the worm and 
directed loop updates [42-46]. In any case, the Monte 
Carlo estimators Eqs. (9,11) are independent to the de¬ 
tailed sampling procedures. It suffices to count kL and 
kR of Monte Carlo configurations to calculate the fidelity 
susceptibility. Examples will be presented in Sec. IV A 
and Sec. IV C. 


2. Stochastic series expansion 

SSE is based on a Taylor expansion of the partition 
function [46], 


B. Implementations 

We now discuss implementation of the estimators 
Eq. (9) and Eq. (11) in various concrete QMC meth¬ 
ods. See Sec. IIB1 for discussions about continuous¬ 
time worldline [42-45] and diagrammatic determinan¬ 
tal [47-53] approaches, and Sec. IIB 2 for discussions 
about stochastic series expansion approach [46]. 


1. Continuous-time worldline and diagrammatic 
determinantal approaches 

Continuous-time worldline methods [42-45] are widely 
used to simulate boson and quantum spin systems, 
while the diagrammatic determinantal approaches are 
the state-of-the-art methods for solving quantum impu¬ 
rity [47-50] and fermion lattice models [51-53]. A com¬ 
mon feature of these methods is to split the Hamiltonian 




n=0 


n\ 



(13) 


which may seem to be different from the framework 
of Eq. (8). However, as is shown in Ref. [56], one 
can formally treat the SSE as the time-dependent ex¬ 
pansion Eq. (12) with respect to the full Hamiltonian 

h = h q + \h 1 . 

In implementation of SSE, one truncates the sum to 
a large number M and pads M — n identity operators 
in the square bracket of Eq. (13). SSE then samples 
operators in the fixed-length operator string. To map 
to a Monte Carlo configuration in the continuous-time 
formalism, one can assign an imaginary-time to each op¬ 
erator as shown in the bottom of Fig. 3. As long as the 
mapping keeps the relative order in the original operator 
string, the Monte Carlo weight remains unchanged [57]. 
In particular, the configuration is sampled with a weight 
proportional to X k if there are k of XH\ operators in the 
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sections. 


A. Ground state 

Using the diagrammatic expansion for the projection 
operator, the unnormalized ground-state wavefunction 
Eq. (10) has the following form, 


Figure 3. Division of the operator string in SSE to measure 
the fidelity susceptibility. The slots represent the fixed-length 
operator string where empty slots hold identity operators, the 
red circles (blue squares) correspond to the operators in Hi 
(Ho). These operators can be mapped to a continuous-time 
configuration indicated by the arrows. A division can then 
be made on the imaginary-time axis, for example at f3/ 2. An 
equivalent approach without explicit mapping to continuous¬ 
time is to divide the operator string at the location indicated 
by the vertical dashed line, where the integer £ is drawn from 
a binomial distribution. For the estimators Eqs. (9,11) one 
counts the number of red circles (operators in XHi) in both 
sides for kh and kR. In this example M — 12, n = 6,£ = 7, 
and kL — kR — 2. 


l*o) 


oo /./3/2 pfi/2 

lim X k / dr \... / dr & x (14) 

J 0 Jr k _i 


(—l) /e e “^/ 2_rfc ^°Hi ... Hie~ Tl ^° 


|®t> • 


Substituting this into Eq. (4), one obtains Eq. (11). 
The estimator also holds for the continuous-time aux¬ 
iliary field expansion methods (CT-AUX [49] and LCT- 
AUX [53]), because one can cast the ground-state wave- 
function to a similar form as Eq. (14), assuming the shift 
parameter used in these methods [51] to be proportional 
to A. 


operator string. In this way, although the sampling of 
SSE is carried out differently from Eq. (12), the general 
framework of Eq. (8) still applies. The fidelity suscepti¬ 
bility is then measured easily by counting the numbers 
kL and kR of operators associated with A Hi in the two 
halves of the imaginary time axis after the mapping. 

From Fig. 3 it is clear that even though one performs 
an equal bipartition in the imaginary-time axis, the cor¬ 
responding location of division is not always in the cen¬ 
ter of the operator string. In fact, it is easier to directly 
sample the location of division in the operator string, as 
shown in the upper part of Fig. 3. A division at the £-th 
position (£ = 0,1,..., M) means that there are i slots 
being mapped to one half of the imaginary-time axis and 
M — i slots to the other half. Therefore the division 
£ itself follows a binomial distribution p(£) = (^Jf) 
which can be sampled directly. In this way, the fidelity 
susceptibility can be efficiently calculated in SSE simi¬ 
lar to the continuous-time QMC approaches discussed in 
Sec. IIBl. As M —> oo, the binomial distribution ap¬ 
proaches to a delta function peaked at the center of the 
operator string. Only in this limit the position in the op¬ 
erator string can be directly interpreted as the imaginary- 
time and a bipartition of the operator string in the center 
will yield the correct result for the fidelity susceptibility. 


III. DERIVATIONS 


B. Non-zero temperature 

We present two derivations of the non-zero tempera¬ 
ture estimator Eq. (9). 


1. Derivation based on the definition of non-zero 
temperature fidelity 

The non-zero temperature estimator Eq. (9) can be ob¬ 
tained directly from the definition of the non-zero tem¬ 
perature fidelity [58], 


F = 


Xr (e-PH(\)/2 e -pH(\+e)/2 


) 


\ (Tr(e-/ 3 *( A ))Tr(e-^( A + £ ))) 


1/2 • 


(15) 


This is a non-zero temperature generalization of Eq. (2) 
and leads to Eq. (7) by using the definition of the fi¬ 
delity susceptibility Eq. (3) [41]. We expand the traces 
of the density matrices around the partition function 
Z = Tr( e -^( A )) to 0(e 2 ), 

Tr (e-^^+O) = Z + ed x Z + jd 2 x Z, 

Tr ^ e -pH(\)/ 2 e - 0 H(\+c)/ 2 ^ = z + e f )xZ + e f-d 2 x Z, 


In this section we derive the estimators for the fidelity 
susceptibility at nonzero temperature Eq. (9) and for 
ground-state projector formalism Eq. (11). Readers may 
skip this section and continue reading with the following 


where the notation d\ indicates that the partial deriva¬ 
tive acts only on operators in the imaginary time interval 
0 < r < 13/2. Substituting the above two expansions into 
Eq. (15) and keeping terms up to 0(e 2 ), one obtains 
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T/O . {BxZf Bjz djz {dxzf 

F 2Z 2 2 Z 4 Z 4 Z 2 

= (fcfl) 2 _ ~ 1)) (fc(fc - 1)) _ (kf 

2X 2 2X 2 + 4X 2 4X 2 

_ {^l^r) — (k L )(k R ) 

2A 2 ' 1 


where k (A) is the number of vertices in the range of 0 < 
r < A. For example, k (/?) = k and k (fij 2) = kR. When 
choosing A = /3 and using G (r) = G (fd — r), Eq. (20) 
becomes 

/3 j\ T G(T) = ^(k(k- 1)). (21) 

When setting A = (3/ 2, Eq. (20) reads 


We have used the partition function in the form of Eq. (8) 
to obtain the second line, and (&l) = (&r) = (k) /2 for 
the third line. This derivation is abstract and is inde¬ 
pendent to the details of a QMC scheme. Carrying out 
a similar procedure starting from Eq. (2), one can also 
prove the ground-state estimator Eq. (11). 


2. Derivation based on the imaginary-time correlator 
Eq. (7) " 

This derivation starts from the definition of fidelity 
susceptibility based on the imaginary-time correlator 
Eq. (7). We utilize its connection to the Monte Carlo 
weight appeared in the Eq. (12) to derive the non-zero 
temperature estimator Eq. (9). First of all, the second 
term in the square bracket of Eq. (7) can be measured 
directly from the average expansion order [46, 47, 50, 51], 


W—W (17) 

Integrating over the imaginary-time and use (ki) = 
(kR) = (k) /2, one has 


rP / 2 
/ dr 

Jo 


~(HiY 


(kf 

8X 2 


(&l) (fcfl) 

2X 2 ' 


(18) 


We then consider the QMC estimator for the first term 
of Eq. (7) 


G(ri— t 2 ) = (T Hi (n) H\ (t 2 ) 


= as ( X S “ Tl ) 5 ( Tj - r 2 ) ) , (19) 




where T is the time ordering operator, r* and Tj are the 
imaginary times of two vertices in the Monte Carlo con¬ 
figuration. Integrating both sides of Eq. (19) and using 
the fact that G (ti — r^) only depends on \t\ — T 2 I, one 
has 



dr 2 G(T 1 - r 2 ) = 


2 / drG (r) (A — r) 

^(fc(A) [fc (A) - 1]), 

( 20 ) 


drG (t) - T^j = ^ {k R {k R - 1)). 


( 22 ) 


Together with Eq. (21) it leads to 


l 




4A 2 
( kjJzR) 
2A 2 * 


2A 2 


(23) 


In combination with Eq. (18) we have derived Eq. (9). 


IV. APPLICATIONS 

We first demonstrate the power of the new approach 
by identifying quantum and classical phase transitions in 
the Bose-Hubbard model and in the spin-1/2 XXZ model. 
Then we use the fidelity susceptibility to address the 
presence of the intermediate quantum spin liquid state 
in the Hubbard model on the honeycomb lattice. In all 
cases this required only minimal modifications to existing 
codes. We have purposely chosen a variety of QMC meth¬ 
ods in the following to demonstrate the wide applicability 
of the covariance estimators. Because of the flexibility of 
the non-zero temperature estimator we used Eq. (9). In 
Sec. V A we compare it to the zero-temperature scheme. 


A. Quantum phase transition in the Bose-Hubbard 
model 

First we use the fidelity susceptibility to probe the 
quantum phase transition in the Bose-Hubbard model, 

- A X (&I&J + b Y b i) ’ 

(id) 

(24) 

where U is the on-site interaction and fi is the chemi¬ 
cal potential. The driving parameter A has the physical 
meaning of a tunneling amplitude. The Bose-Hubbard 
model has a well-known quantum phase transition be¬ 
tween the Mott insulating state and the superfluid state 
as X/U increases [1, 2]. In particular, for integer fill¬ 
ings the system has an emergent Lorentz invariance at 
the critical point and the dynamical critical exponent is 
z = l[7]. 


H ~ X ( - 1) - 
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The fidelity susceptibility has previously been calcu¬ 
lated using density-matrix-renormalization-group for the 
one-dimensional Bose-Hubbard model [60-62]. We now 
calculate the fidelity susceptibility on a square lattice 
with N = L 2 sites at unit filling by tuning fi. In ac¬ 
cordance with the dynamical critical exponent z = 1, 
we scale the inverse temperature proportionally to the 
system length (3U — L. The simulation employs the di¬ 
rected worm algorithm [43, 63, 64]. We utilize Eq. (9) to 
sample the fidelity susceptibility by counting the num¬ 
ber of kinks in the worldline configuration, as illustrated 
in Fig. 2(a). Figure 4 shows that as the system size 
increases, the peak in the fidelity susceptibility (as a 
function of the driving parameter A) is becoming more 
pronounced and is shifting towards the previously de¬ 
termined critical point (A /U) c = 0.05974(3) [59]. The 
ability to calculate the fidelity susceptibility using the 
state-of-the-art directed worm algorithm [43, 63, 64] will 
greatly advance the study of quantum phase transitions 
of ultracold bosons. It is worth to point out the fidelity 
susceptibility is related to the quantity (kinetic-energy 
correlator) previously calculated in the study of Higgs 
mode in a two-dimensional superfluid [65]. 


B. Classical phase transition in the XXZ model 

Next we consider the spin-1/2 antiferromagnetic XXZ 
model on a square lattice with N = L 2 sites, 

A = J *H + A £ (sf S* + S?§y ) , (25) 

(hi) (i,o) 

where the driving parameter A plays the role of the cou¬ 
pling strength in the XY-plane. When A dominates the 
Hamiltonian favors Neel order in the XY-plane, while if 



Figure 4. Fidelity susceptibility per site of a Bose-Hubbard 
model on a square lattice at unit filling. The vertical line 
indicates the critical point determined in Ref. [59]. 


J z dominates the system has an antiferromagnetic Ising 
ground state. The Heisenberg point A = J z is a quan¬ 
tum critical point, which separates the XY order and the 
Ising order. This quantum critical point can be easily 
located from the peak of the fidelity susceptibility (not 
shown). Our approach makes it possible to obtain the fi¬ 
delity susceptibility in much larger systems compared to 
the previous exact diagonalization study [67], thus can 
enable a more accurate scaling analysis. 

At nonzero temperature, thermal fluctuations will de¬ 
stroy the antiferromagnetic Ising phase at a second-order 
phase transition. Since one can cross the phase bound¬ 
ary either by changing A or the temperature T, we see 
that the fidelity susceptibility can also indicate thermal 
phase transitions. As a demonstration we fix A = 1, 
J z = 1.5 and scan the temperature T to drive a phase 
transition from the low-temperature antiferromagnetic 
Ising phase to the high-temperature disordered phase. 
Figure 5 shows the fidelity susceptibility calculated using 
Eq. (9) via the SSE method [46, 68]. The peak in the 
fidelity susceptibility correctly single out the previously 
determined critical temperature (T/ A) c « 0.75 [66]. 

C. Intermediate phase in the Hubbard model on 
the honeycomb lattice 

Finally we apply the fidelity susceptibility estimator 
to a more challenging and controversial example - the 
Hubbard model on the honeycomb lattice, 

H = Y (4Av + C^Ciaj 

(i,o) o‘={t,4} 

+ A Y ( r V - 0 («4 - 0 > ( 26 ) 



Figure 5. Fidelity susceptibility per site of a XXZ model on 
square lattice versus temperature. The vertical line indicates 
the critical temperature determined in Ref. [66]. 



















Figure 6. Fidelity susceptibility per site of the Hubbard model 
on the honeycomb lattice Eq. (26) with N = 2 L 2 sites. 


where A has the meaning of on-site Hubbard interaction 
strength. The simulation employs the recently devel¬ 
oped efficient continuous-time QMC method for lattice 
fermions (LCT-INT) [53] [69]. We consider lattices with 
N = 2 L 2 sites, with L = 6, 9,12 and scale the inverse 
temperature fit = L. 

The ground-state phase diagram of the Hubbard model 
on the honeycomb lattice [70] has been controversial. It 
was suggested to possess an intermediate non-magnetic 
spin-liquid phase for A jt E [3.5,4.3] [71]. However, more 
recent QMC studies on larger systems [72] and with im¬ 
proved observables [73, 74] suggest a single continuous 
phase transition at X/t ~ 3.8 belonging to the Gross- 
Neveu universality class [75]. Other less unbiased meth¬ 
ods such as quantum cluster approaches give conflicting 
results on the presence of the intermediate phase [76-81], 
depending on implementation details. 

The fidelity susceptibility offers a new perspective on 
the debate about the phase diagram. In the scenario 
with an intermediate phase, there shall be two features 
in xf when X/t approaches the two phase boundaries. 
This consideration is independent of the presence of a 
local-order-parameter description of the possible inter¬ 
mediate phase. Figure 6 shows the fidelity susceptibility 
per site for various system sizes obtained using Eq. (9). 
It exhibits a single broad peak for small systems (and 
high temperature). The peak becomes sharper and shifts 
towards smaller interaction strength as the system size 
increases. The fidelity susceptibility data presented in 
Fig. 6 is consistent with a single phase transition at 
X/t ~ 3.8. In future studies for larger system sizes, 
and also with possible extension to a continuous range 
of A (by using histogram reweighting [82, 83] or quan¬ 
tum Wang-Landau approaches [84]), it may be possible 
to precisely determine the critical point, or even the crit¬ 
ical exponent solely from the fidelity susceptibility data. 
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r/P 


Figure 7. (a) Histogram obtained by counting separable ver¬ 
tex pairs with distance r compared with the exact result of 
(Hi (r) Fi)xmax{r, ft— r}. (b) Histogram obtained by count¬ 
ing vertices with distance r compared with the exact results of 

(Hi (t) Hi). For r close to /3/2, the correlator approaches to 

~ 2 

(Hi) , indicated by the dashed blue line. These simulations 
are performed for the Hubbard model Eq. (26) on a four-site 
open chain with X/t = —2 and fit = 8. 

V. DISCUSSION 

To help the reader gain a better understanding of the 
estimators Eq. (9) and Eq. (11), we first discuss their 
relationship then compare them with the previous ap¬ 
proach adopted in SSE calculations [24, 25]. Finally, we 
compare the fidelity susceptibility approach with other 
generic approaches for detecting phase transitions. 

A. Relation of the ground-state and non-zero 
temperature estimators 

The factor of two difference in Eq. (9) and Eq. (11) 
is due to the different boundary conditions of the 
imaginary-time axis in the ground-state projection and 
non-zero temperature QMC formalisms, see Fig. 1. We 
use a four-site Hubbard model (Eq. (26)) as an illustra¬ 
tive example. Consider the integrand of Eq. (7), the 
correlator G(r) = (Hi(r)Hi) is related to the distri¬ 
bution of the vertices on the imaginary-time axis. For 
a given configuration, the probability of finding two ver¬ 
tices with a time difference r is proportional to A 2 G(r). 
If we equally divide the imaginary-time axis into two 
halves and impose the additional constraint that the 
two vertices reside in different halves (denoted as a 
separable vertex pair), the joint probability changes to 
A 2 G(r) min{r, ft — r}. Figure 7(a) shows the histogram 
of separable vertex pairs accumulated in the imaginary 
time, which indeed agrees with the exact curve. Sum¬ 
ming up the histogram gives the total number of separa- 
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Figure 8. QMC results for the fidelity susceptibility com¬ 
pared with exact results (solid line). The non-zero temper¬ 
ature QMC data (red dots) is obtained from Eq. (9), while 
the ground-state data (blue square at (3~ 1 — 0) is obtained 
from Eq. (11) in a projector LCT-INT calculation [85]. The 
system is the same as in Fig. (7). 


ble vertex pairs, which equals to the following integration, 
rP /2 rP 

(&l&r) = A 2 / drG(r)r + A 2 / drG(r)(f3 — r). 

Jo J(3/2 

(27) 

Since G(r) is symmetric around r = /3/2 in the non¬ 
zero temperature simulation, the two terms of Eq. (27) 
are equal. Thus Eq. (27) reduces to Eq. (23). Further¬ 
more, Fig. 7(b) shows the correlator G(r) sampled by 
accumulating the histograms of distances between ver¬ 
tices [48, 50] together with the exact results (solid black 
line). The correlation between vertices decays rapidly 

with imaginary-time distance and approaches to the un- 

^ 2 

correlated value (Hi) (dashed blue line). 

However, in the zero-temperature limit, the correlator 
G(r) decays monotonically with r and two vertices will 
decorrelate for r > /3/2, where /? —» oo in the projection 
scheme. Therefore the second term of Eq. (27) reduces to 

A ^ and cancels half of the second term in 

the estimator Eq. (11), resolving the apparent difference 
by the factor of two. In practical calculations, it is how¬ 
ever crucial to adopt the correct formula to obtain consis¬ 
tent results, as illustrated in Fig. 8. The fidelity suscepti¬ 
bility calculated using Eq. (9) in a non-zero temperature 
LCT-INT [53] simulation agrees perfectly with exact di- 
agonalization results. The blue square shows the value 
obtained using Eq. (11) in a projector LCT-INT calcu¬ 
lation [85], which correctly reproduces the exact value of 
the ground-state fidelity susceptibility. 

Figure 7(b) also reveals the difficulty of computing the 
fidelity susceptibility. If the decay of G(r) is faster than 
1/t, the integrand of Eq. (7) has vanishing contributions 
at large r. However, as the two terms in Eq. (7) are sam¬ 
pled independently in the actual QMC simulations, un¬ 
correlated vertices at large imaginary-time distance will 
cause noises in the fidelity susceptibility signal. For the 
applications in Sec. IV we thus perform the calculations 
at nonzero temperature as it provides a natural cutoff. 


Figure 9. The weight function according to Eq. (29) for vari¬ 
ous truncation lengths M in a SSE calculation. 


B. Comparison to previous approaches 


The present approach to sample the fidelity suscepti¬ 
bility is more generic and efficient than those developed 
in Refs. [24, 25] specifically for the SSE method. It is nev¬ 
ertheless instructive to compare them in detail. The key 
difference lies in the sampling of the first term of Eq. (7). 
References [24, 25] employ the SSE estimator [87] [88] 


G(r) 


M - 1 
A 2 /3 2 




(<?(«)>, 


where G(n) is the number of occurrences of two operators 
from Hi that are separated by n positions in the fixed- 
length operator string (n = 0 if they are next to each 
other). Multiplying both sides with max{r, (3 — r} and 
integrating over the imaginary time, one finds 


rP 

A 2 / dr G(t ) max{r, p 

Jo 


M—2 


~ r} = W(n) ( G(n )), 


n =0 


(28) 


where the weight function W(n) is written in terms of 
the regularized incomplete beta-function I x (a,b) [89], 


77 T 1 

W(n) = Ii(n + 2,M-n - 1)—V- + 

' 2 v M 

Ii(M~n,n + 1) — ^ (29) 

References [24, 25] explicitly go through k(k — l)/2 
pairs of vertices to accumulate ( G(n)) and multiply it 
with the weight function W(n). In Eq. (27), however, 
the multiplication by the imaginary time r is taken into 
account implicitly by the sampling procedure (which re¬ 
quires separable vertices). Besides being more generic, 
our approach reduces the computational cost from 0(k 2 ) 
to 0(k), which is crucial for the simulation of bosonic and 
quantum spin systems. In this sense, the specification of 
our general result Eq. (9) for the SSE method can be 
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regarded as an improved estimator of Eq. (28), which 
by itself already improves the approach of Refs. [24, 25] 
in several aspects [90]. The improved estimator Eq. (9) 
not only unifies the SSE approach in a broader context 
of continuous-time diagrammatic QMC methods, it also 
gives better statistics with less computational cost com¬ 
pared to Eq. (28). 

Figure 9 shows the weight function W(n) for various 
truncation lengths. As M increases, it approaches to two 
straight lines, and a division in the center of the opera¬ 
tor string would yield increasingly accurate result for the 
fidelity susceptibility, consistent with the discussion in 
Sec. IIB 2 concerning the large M limit. 


C. Relationship to other quantities 

The fidelity susceptibility is related to the second-order 
derivative of the free energy A = — ^ In Z [25,91]. Due to 

the Helhiiann-Feynman theorem [92, 93], (Hi) equals to 
the first-order derivative of the free energy with respect 
to A. A further derivative following the Kubo formula 
gives, 


d 2 A _ <9 (Hi) 
~d)J ~ d\ 



dr 


(Hi (t) Hi) — (HiY 


Of) z Off - Of 

-/?A 2 


(30) 


The third equality follows from Eq. (17) and Eq. (21) [94]. 
The quantity resembles the widely used SSE estimator 
for the specific heat [46, 95], but can be used to probe 
quantum phase transitions [25]. At zero temperature 

= ~ X d ^d\* ’ an d latter quantity was computed 
using numerical differentiation of the kinetic energy, so 
as to address the quantum phase transition in the Hub¬ 
bard model on the honeycomb lattice [71, 81]. As is 
pointed out in Refs. [25, 91], the fidelity susceptibility 
has a stronger singularity compared to the second-order 
derivative of the free energy and is thus a better indicator 
of quantum phase transitions. There are concrete exam¬ 
ples in a class of topological phase transitions, which do 
not exhibit singularity in the second-order derivative of 
the ground-state energy [96], but can still be detected 
using the fidelity susceptibility [97]. 

The covariance which appeared in the estimators 
Eq. (9) and Eq. (11) can also be written as 

(k L k R ) - (k L ) (, k R ) = i [Var(fc) - Var(fc L ) - Var (k R )], 

where Var(x) = ( x 2 ) — ( x ) 2 is the variance of x. This 
expression has an appealing meaning, i.e., the distribu¬ 
tions of the vertices residing on the whole and on the 
halves of the imaginary-time axis have different widths, 
and the difference in these widths gives the estimate. 
This form resembles the bipartite fluctuation [98], which 


was proposed to be a diagnostic tool for phase transi¬ 
tions [99] because of its relation to the entanglement en¬ 
tropy. However, there are important differences. First, 
the fidelity susceptibility estimator requires a division in 
the imaginary-time axis for vertices, not in the real space 
for the physical particles. Second, the total number of 
vertices is fluctuating in the QMC simulations as opposed 
to being conserved in the case of bipartition fluctuations. 
Third, it is easier to locate the critical point using the 
fidelity susceptibility. As is shown in this paper and in 
many previous studies [16], the fidelity susceptibility ex¬ 
hibits an increasingly sharp peak at a phase transition as 
the system size enlarges. On the other hand, to utilize bi¬ 
partite fluctuations and entanglement entropy for phase 
transition, one typically needs to resolve the scaling or 
subleading behavior with the system size, which is often 
difficult in finite size simulations. 


VI. OUTLOOK 

We have presented a general approach to compute the 
fidelity susceptibility of correlated fermions, bosons, and 
quantum spin systems in a broad class of quantum Monte 
Carlo methods [42-53]. The calculation of the fidelity 
susceptibility is surprisingly simple yet generic. It pro¬ 
vides a general purpose indicator of quantum phase tran¬ 
sitions without the need for a prior knowledge of the local 
order parameter. 

Conceptually, our work shows it is rewarding to view 
the modern QMC methods [42-53] in a unified frame¬ 
work provided by Eq. (8) , which deals with the same type 
of classical statistical problem irrespective of microscopic 
details of the original quantum system. In the QMC sim¬ 
ulations, a quantum phase transition manifests itself as 
a particle condensation transition driven by changing of 
the fugacity of the corresponding classical model. This 
connection suggests generic ways to detect and charac¬ 
terize quantum phase transition through studying classi¬ 
cal particle condensations. For example, Eq. (30) actu¬ 
ally relates the second-order derivative of free energy of a 
quantum system to the particle compressibility of a vir¬ 
tual classical system. In this respect, the significance of 
the covariance estimators Eqs. (9,11) is evident because 
they capture the key critical fluctuation upon a particle 
condensation transition. 

It is straightforward to generalize our Eqs. (9,11) to 
cases with multiple driving parameters, where one needs 
to count the vertices of different types (as is already done 
in the SSE calculations in Sec. IIB 2 and Sec. IV B). It 
is interesting to find out whether this can lead to a gen¬ 
eral approach to measure the Berry curvature (the imag¬ 
inary part of the quantum geometric tensor) in quan¬ 
tum Monte Carlo simulations. Related to these efforts, 
the non-equilibrium QMC method is developed in recent 
years [100, 101] to study non-adiabatic response of quan¬ 
tum systems in the imaginary time. In particular, it also 
allows the extraction of the fidelity susceptibility and the 
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Berry curvature [100, 102]. It would be interesting to 
compare the non-equilibrium QMC approach [100, 101] 
to the equilibrium one presented in this paper. 

Last but not least, the Hamiltonian Eq. (1) has further 
implications beyond quantum phase transitions. The ef¬ 
ficient estimators Eqs. (9, 11) may also provide useful 
insights in the simulations of adiabatic quantum com¬ 
putation [103-105] and non-adiabatic quantum dynam¬ 
ics [106, 107]. 
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